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ABSTRACT 

Cooling and heating functions of cosmic gas are a crucial ingredient for any study of gas dynamics and 
thermodynamics in the interstellar and intergalactic medium. As such, they have been studied extensively in 
the past under the assumption of collisional ionization equilibrium. However, for a wide range of applications, 
the local radiation field introduces a non-negligible, often dominant, modification to the cooling and heating 
functions. In the most general case, these modifications cannot be described in simple terms, and would require 
a detailed calculation with a large set of chemical species using a radiative transfer code (the well-known code 
Cloudy, for example). We show, however, that for a sufficiently general variation in the spectral shape and 
intensity of the incident radiation field, the cooling and heating functions can be approximated as depending 
only on several photoionization rates, which can be thought of as representative samples of the overall radiation 
field. This dependence is easy to tabulate and implement in cosmological or galactic-scale simulations, thus 
economically accounting for an important but rarely-included factor in the evolution of cosmic gas. We also 
show a few examples where the radiation environment has a large effect, the most spectacular of which is a 
quasar that suppresses gas cooling in its host halo without any mechanical or non-radiative thermal feedback. 

Subject headings: methods: numerical 



1. INTRODUCTION 

The ability of cosmic gas to radiate its internal energy (i.e. 
radiative cooling) and to absorb energy from the incident radi- 
ation field (radiative heating) is a primary distinction between 
the gas and dark matter; radiative heating and cooling pro- 
cesses are important in almost every study of gas dynamics 
or thermodynamics in the interstellar and intergalactic media. 
Because of this importance, cooling processes in the gas have 
been investigated in numerous prior studies, appear as central 
chapters in multiple textbooks, and are computed by several 
publicly available codes. 

However, while the physics of radiative cooling and 
heating is well understood, the actual application of cooling 
and heating functions for studies of interstellar and inter- 
galactic gas is surprisin gly incomplete. The classic "standard 
coolin g function" (e.g. ICox & Tucker! 1 19691: iRavmond et all 
1976b IShull & van Steenberd 119821; iGaetz & Salpeterl 
1983[ iBoehringer & Hensleri Il989t ISutherland & Dopit 
1991 




ISht 

lLandi & Landinil Il999t 



Benjamin e t alj 12001 



Santoro & Shulll 120061: iGnat & Sternberg! 120071: ISmith et alJ 
2008TT has indeed been computed and tabulated quite pre 



cisely. However, the "standard cooling function" is computed 
under the assumption of pure collisional ionization equi- 
librium (CIE), which is not always valid in the interstellar 
medium and is never valid in the intergalactic medium (c.f. 
iWiersma et al.l 120091) . In many astrophysical applications 
the incident radiation field introduces significant, often 
dominant, modifications to the "standard cooling function". 
On top of that, in some environments the assumption of the 
photoionization equilibrium may not be sufficientl y accurate 
(ISutherland & Dopital!T99l I Santoro & Shull2006l) . 
Such dependence can be illustrated by comparing the pure 
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FIG. 1 . — Example of the importance of the incident radiation field on the 
cooling functions: blue dashed and solid lines show the "standard cooling 
function" for the metal-free and solar-metallicity gas respectively. Corre- 
sponding red lines show the same cooling functions for the fully ionized gas. 

CIE cooling function with the cooling function in fully ion- 
ized gas, as shown in Figure Q] for both metal-free and solar- 
metallicitjQ gas. In the fully-ionized limit, where the only 
cooling process is bremsstrahlung, the cooling function over 
a wide range of temperatures differs from a pure CIE case by 
more than two orders of magnitudel 

Thus, the cooling function for photoionized gas depends not 
only on the gas temperature, number density, and metallicity, 
but also on the incident radiation field. There is, of course, 
nothing new in that statement. The crucial role of the radia- 
tion environment has always been understood by practition- 
ers in the field. The challenge, however, is in economically 
accounting for this dependence in full 3D numerical simula- 
tions, where the cooling and heating functions are evaluated 
billions or even trillions of times during a single simulation. 
In the worst case scenario - the radiation field /„ varies ar- 

4 Throughout this paper, "solar metallicity" refers to the metallicity of the 
gas in the solar neighborhood, Z ss 0.02, not the actual metallicity of the Sun. 
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bitrarily - one can introduce sharp edges and features in the 
radiation spectrum that are specially designed to ionize par- 
ticular levels of particular elements. This allows the cooling 
function to be "sculpted" in an essentially arbitrary way. 

One possible way to account for the effect of the inci- 
dent radiation field is to fix the radiation spectrum and am- 
plitude. For example, in studies of intergalactic medium it 
is often (but not always) sufficient to account for the cosmic 
background radiation. Since the cosmic background radia- 
tion evolves with redshift, the cooling and heating functions 
become redshift-dependent, but such 1 -dimensional depen- 
dence is easy to pre-compu t e and tabulate f or use in sim- 
ulations dBensonet all 120021: iKravtsovllIOOll: IWiersma etafl 
120091: 1 Vasilievll20 1 ll) . Unfortunately, these cooling and heat- 
ing rates are then often used for modeling gas dynamics in 
galactic halos or even ISM - environments where the cosmic 
background radiation is a sub-dominant component of the in- 
cident radiation field. 

Therefore, it is desirable to find a way to account for a gen- 
eral shape of the incident radiation field without the need to 
recompute the cooling and heating functions every time they 
are needed. In this paper, we show that it is possible to come 
up with an approximate solution for this problem using a suf- 
ficiently general model for the radiation field spectrum. 

2. APPROXIMATING THE COOLING AND HEATING FUNCTIONS 

The radiative term in the internal energy equation - the rate 
of change of the gas internal energy due to radiative gains and 
losses - can be represented as 



dU 
dt 



= n 2 b [T(T,...)-A(T, ...)], 



(1) 



lad 



where U is the gas thermal energy and n b = rin + 4rin e + .. is 
the total baryon number density. We explicitly factored out n\ 
in both the cooling (A) and the heating (T) functions so that 
these are density-independent in the CIE limit. 

In the most general case the cooling and heating functions 
depends on an extremely large set of arguments: gas temper- 
ature T, baryon number density n b (in addition to rr b depen- 
dence explicitly accounted for in Equation ([T]i), the fractional 
abundance Xy for the species i (including atomic and ionic 
species, various molecules, and cosmic dust) at level j, the 
distribution of the column density for the species i at level j 
at different velocity values with respect to the systemic veloc- 
ity dNjj(v)/dv, the specific intensity of the radiation field as a 
function of frequency /„, and the heating rate by cosmic rays 

CcRi 



J-(l , ...) - J- I l ,n h ,Ajj, — — — ,./„,C,cr 
where hereafter T denotes either T or A, 



(2) 



r(...) 

AC) 



Obviously, such a complex dependence cannot be described 
in simple terms, and would require a detailed calculation with 
a large set of chemical species using a radiative transfer code 
- for example, the well-known code Cloudy dFerland et al.l 
1 1 9981) . That would make it impractical as a method for com- 
puting the cooling and heating function in realistic three- 
dimensional numerical simulations. 

We, therefore, adopt several major simplifications. First, we 
restrict our focus to a purely optically thin case (all N,j = 0). 



Second, we exclude cooling and heating due to molecules, 
dust, and cosmic rays, since these processes crucially depend 
on radiative transfer and computing them in the optically thin 
limit does not make much physical sense. 
With these restrictions, Equation becomes 

HT,...) = J r (T,n b ,X ij ,J v ). 

Even this is way too complex, as the cooling and heating func- 
tions depend on hundreds of individual level populations for 
atomic and ionized species. 

In the next simplification step we assume that all atoms and 
ions are in the ionization equilibrium, and the level popula- 
tion is in the equilibrium as well. This assumption is actu- 
ally valid in a vast majority of astrophysical environments. In 
the limit of ionization equilibrium, the cooling and heating 
functions only depend on the total abundance of each chem- 
ical element. Finally, if we assume that the abundance pat- 
tern for heavy elements is solar, and ignore small variations 
of helium abundance, then the cooling and heating functions 
become just functions of the gas metallicity, 



F(T,...) = F(T,n b ,Z,J v ). 



(3) 



Often, Equation ([3]) is what is actually called "cooling" and 
"heating functions". For example, the CIE cooling and heat- 
ing functions are just 

F cm (T,Z) = F(T,n b ,Z,0) 

(which, in this limit, is also independent of n b ). 

At low enough densities and faint enough incident radiation 
fields, most of reactions that result in cooling and heating in 
gas are interactions of an atom/ion with either a photon or an 
electron. Hence, in this limit cooling and heating functions 
(Eq.Q} can be substantially simplified: 



F(T,n b ,Z,J v )K, T{T,Z,—) 
n b 



Unfortunately, this approximation is only valid for rather 
low values of the radiation field; for example, it is only 
marginally valid in typical ISM conditions in the Milky Way. 
Several physical processes break the ideal density-squared 
dependence. At high enough densities various 3-body pro- 
cesses become important - in particular, 3-body recombina- 
tion can become important at low temperatures for densities 
as low as 10 -4 cm" 3 . For hard enough incident radiation spec- 
tra secondary ionizations and heat deposition from secondary 
electrons introduce complex density dependence. For strong 
enough radiation fields some of highly excited energy levels 
have critical densities within the density range we consider 
here (n b < 10 6 cm~ 3 ). 

All of these processes, however, are relatively smooth func- 
tions of the density at a constant value of J u /n b . Hence, with- 
out any loss of generality, we can re-write Equation (0) as 



F(T,...)=F(T,Z,— ,n b ). 

n b 



(4) 



The advantage of this transformation is that the third argument 
includes most of the density dependence; the explicit density 
dependence parametrized by the fourth argument is relatively 
weak and can be accounted for in a numerical implementation 
in an economic manner (the fact we take a full advantage of 
below). 
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Finally, we exclude from our cooling and heating functions 
Compton cooling and heating and free-free absorption - not 
because they violate the ansatz (O, but because we found 
empirically that the functional dependence of those two pro- 
cesses on the properties of the medium is sufficiently differ- 
ent from the characteristic dependence of line excitation and 
emission, so that numerical approximations that we discuss 
below become substantially less accurate with those two pro- 
cesses included. In addition, Compton processes and free-free 
absorption depend on the whole shape of the radiation spec- 
trum, even in the low energy regime that is unimportant for 
the line excitation and ionization. These two processes are de- 
scribed by sufficiently simple analytical expressions that can 
be easily added separately to a numerical code if there is such 
a need. 

3. MODELING COOLING AND HEATING FUNCTIONS AND THE 
INCIDENT RADIATION FIELD 

As we mentioned above, it is not possible to account para- 
metrically for an arbitrary radiation field spectrum. However, 
in many astrophysical applications the incident radiation field 
is dominated by radiation from stars, AGN or a combination 
thereof. Thus, we model the incident radiation field as 



J v = Joe T " 



1 



1+/q- ■ l+/ G 



-Su + 



(5) 



where x is the photon energy in Rydbergs (x = hv/{ \ Ry)), 
and s„ is a fit to t he stellar spectrum from Starburst99 
dLeifherer et al.|[l999l) . 



1 



'5.5, 

-ii 



5.5 1 0.4x" LI 



x < 1 

1 <x< 2.5 
2.5 <x <4 



2x 10-V/(exp(x/1.4)-l), A<x 



(this fit is shown in Fig. 4 of iRicotti et al.l (120021) '). Equation 
([5]) also includes the possibility that the incident radiation field 
is attenuated by gas with the opacity 



CHI.0 



[0.76osi(f)+0.06osei(f)] J 



where <Jj{v) and cr^o are photoionization cross-sections and 
their values at respective ionization thresholds for j = HI and 
He I, and tq is a parameter. 

Overall, the radiation field model from Equation (O con- 
tains 4 parameters: the amplitude Jo, the AGN-like power- 
law contribution slope a, the ratio of the AGN-like to stel- 
lar component /g, and the shielding optical depth parameter 
To- The last 3 parameters are dimensionless; we choose to 
measure Jo in units of the typical radiation field in the Milky 
Way galaxy, Jmw = 10 6 p hotons cm" 2 s" 1 ster" 1 eV -1 dDraind 
[r^lMathisetaD[l983h . 

For each set of paramet ers, we use the wid ely known pho- 
toionization code Cloudy (Fer land et al.ll 19981) to compute the 
cooling and heating function for a range of gas temperatures 
at fixed gas density and metallicity. Examples of such com- 
putations are shown in Figure |2] For all 3 cases the radiation 
field is the same at 1 Ry, but differs in spectral shape at other 
frequencies (a stellar spectrum, a power-law spectrum, and a 
power-law spectrum shielded by a tq = 100 cloud). In order to 
enforce the optically thin case, we restrict Cloudy calculations 
to a single zone of a negligibly small size. 



We explore the cooling and heating functions for our radia- 
tion field model by sampling the full parameter space (metal- 
licity, density, and the radiation field) on the following grid of 
values: 

Z/Z Q = 0,0.1,0.3,1,3 



lg(n b / cm 3 ) = -6,-5,...,6 
lg(/ cm" 3 jn h /J M w) = -5 , -4 .5 , -4 , . . . , 7 

a = 0,0.5, 1,1.5,2,2.5,3 
lg(/,) = -3,-2.5 -2,...,1 
lg(r ) = -l -0.5,0,..., 3 



(6) 



This parameter range is wide enough to include both extremes 
shown in Fig. Q] the case where the radiation field is com- 
pletely negligible and the case where the gas is fully photoion- 
ized. 

For each of the 5 x 13 x 25 x 7 x 9 x 9 920,000 sets of 
parameters from this grid, we run Cloudy to compute the cool- 
ing and heating functions for 8 1 values of the temperature be- 
tween 10 K and 10 9 K in steps of 0.1 dex (almost 75 million 
Cloudy runs altogether). Using this large database, we now 
consider the various dependencies of the cooling and heating 
functions one by one. All of our subsequent approximations 
are extensively tested below in §0] 

3.1. Metallicity and Density Dependence 

In a further simplification, we expand both cooling and 
heating functions into Taylor series in metallicity up to the 
quadratic term, 



Z 

Z^ 



?2, 



(7) 



where all functions J 7 , depend only on T, J^/n/,, and n b . 

We achieve this decomposition in practice by fitting a sec- 
ond degree polynomial to the five Z values that we sample 
in Table ©. The error introduced by dropping cubic and 
higher power terms is by far the smallest of the errors intro- 
duced by our approximations - in the rms sense, the second 
order expansion of the Taylor series is accurate to better than 
3% - as long as we restrict Z to less than 3 solar metallici- 
ties. The quadratic approximation rapidly loses accuracy as 
the metallicity increases. At metallicities above 5Z Q , approx- 
imation (|7]i even results in negative cooling functions in a few 
instances. 

Six functions T, and A, (i = 0, 1,2) can be used directly, 
but since cooling and heating functions are not necessarily 
monotonic functions of Z, some of T\ and Ti (again, T stands 
for either T or A) can be negative. Since interpolation in log- 
log space is usually more accurate than direct interpolation, 
positive functions are much more suitable for tabulation and 
interpolation. Hence, we replace 6 functions T; with 6 new 
functions J 7 ,- as 

•F\ = •> o + -F\ + -Fit 
T 2 = Fo + 2Fx+AT 2 , 

where symbol T also means either the cooling or the heating 
function. Functions J 7 , are none other than the cooling and 
heating functions at Z = i x Z Q and hence are always positive. 
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FIG. 2. — Incident radiation fields (top panels) and gas cooling (blue lines) and heating (red lines) functions (bottom panels) for three different radiation field 
models: [J ,a,f Q ,T ] = [10O/ M w,3, 10 3 ,0.1] (left), [10O/ M w,2, 10,0.1] (middle), and [1007 M w,2, 10, 100] (right), all forn H = 1cm" 3 . As inFig.Q] solid and 
dashed lines are for Z = Zq and Z = respectively. It is clear that the cooling and heating functions are strongly dependent on the incident radiation field. 



The transformation between J 7 ,- and Ti is linear and can be 
trivially inverted. 

In the following, we always operate on functions J 7 , and 
convert them back to J 7 ,- (i.e. T, and A,) as the very last step. 

3.2. Radiation Field Dependence 
So far we still have not resolved the main challenge - the 
fact that the 6 functions T; that we need to describe depend 
on the whole incident radiation field /„, 

fi = Tj(T,—,n h ). 
n b 

The primary contribution of this paper is that we further ap- 
proximate this dependence by replacing the full radiation field 
with a finite set of photoionization rates. 

Specifically, let us define a normalized rate Qj as 



Qj 



rib 



where Pj is a photoionization rate for some atom or ion, 

Gj(y)n v dv, 



Pj = C 



where cry is the photoionization cross-section and n v is the 
radiation field expressed as the number density of photons at 



the frequency v. We now seek an approximation of the form 



Jv 

Fi{T,—,n b ) « Fi(T,Qj,n h ) 
n b 



(8) 



for i = 0, 1,2 and some set of Qj. If such an approximation is 
possible, then a user would only need to compute the (normal- 
ized) photoionization rates Qj from his/her assumed radiation 
field (provided, such a field is close enough to our assumed 
shape - as we mentioned above, it is always possible to sculpt 
the cooling and heating functions by appropriate choosing the 
radiation field spectrum). 

We present our specific choice for Qj in ■ 

3.3. Explicit Density Dependence 

Finally, we need to address the remaining density depen- 
dence in Equation ©. The trick of using Qj makes this depen- 
dence relatively weak, although highly non-trivial. We adopt 
the simplest possible approach - we tabulate J 7 ,- at the 13 den- 
sity values we tested in Table © and linearly interpolate in 
log-log space. The guaranteed positiveness of J 7 , becomes 
crucial when working in logarithmic space. 

We verified that the error introduced by such interpolation 
is completely sub-dominant to the error of our main approxi- 
mation (O. 

3.4. Notes on the Specific Implementation 
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It makes sense that the rates we choose to represent the 
radiation field should sample the wide range of frequencies. 
For example, since CII is an important coolant in the low- 
temperature regime, one of the rates should sample the ra- 
diation field below the hydrogen ionization threshold. We 
choose the photodissociation rate of molecular hydrogen in 
the Lyman-Werner band as such a rate, simply because that 
rate is also useful for several other processes that can be mod- 
eled in the numerical code (for example, the destruction of 
molecular hydrogen). It also makes sense to use the hydrogen 
and helium ionization rates since both elements are impor- 
tant coolants at T > 10,000 K for all but the highest radiation 
fields. Finally, one of the selected rates should be sensitive to 
high energy photons. 

While it is not possible to explore the full set of some 
600+ photoionization rates for common chemical elements, 
we searched the full hydrogen-like sequence all the way to 
FeXXVI for the most accurate approximation for the cooling 
and heating functions. 

For a practical numerical implementation, a table with cool- 
ing and heating functions should not exceed about 1 GB of 
memory, and a much smaller memory footprint is much pre- 
ferred. Most modern supercomputers offer between 1 and 2 
GB of memory per computing core; a 1 GB table would leave 
no memory for other data structures in a pure MPI code on a 
1 GB/core machine. While pure MPI codes are getting more 
and more rare, even a hybrid code that uses just one MPI task 
per an 8-core node with 8 GB of memory would encounter 
difficulty using a table in excess of 1 GB. 

With these constraints in mind, we explored a wide range 
of possible 3- and 4-dimensional tables. We have not found a 
4-dimensional solution that is sufficiently superior to our best 
3-dimensional table and that fits within our imposed 1 GB 
memory limit. 

Hence, as our primary implementation, we present a 3- 
dimensional table that is constructed from four normalized 
photoionization rates, 2lw, Q m , Qu e i, and Qcvi, combined 
into 3 independent parameters, 




0.353 



-0.103 



0.923 



-0.375 



9l 

Ql 



CVI 



0.263 



Ql 



Ql 



cvi 



0.976 



(9) 



There exist several other combinations of various rates that 
result in almost equivalent approximations. For example, 
adding Qivigxn as a fifth rate reduces the error of the approx- 
imation by about 0.03 dex - not significant enough, in our 
opinion, to justify computing an extra photoionization rate. 

We implement the approximation (O by constructing a grid 
of rj values and computing the average cooling and heating 
functions for all incident radiation fields that happen to have 
the same values for rj. Specifically, we use a logarithmically- 
spaced table for-14.5 < lg(n) < -3, -9.5 < lg(r 2 ) < 0.5, and 
~8 < lg(/3) < ~0-5 with the logarithmic step of 0.5 dex. We 
found that using a finer step in the table does not lead to any 
increase of accuracy. Such a table includes 24x21 x 16 = 
8064 entries; each entry contains a sub-grid of 13 density val- 
ues by 8 1 temperature values, with 6 numbers J 7 , at each grid 
point. A full table takes about 192 MB of memory. The table 
can be further compressed by eliminating entries with similar 
cooling and heating functions, to the total of 3795 entries (92 




T (K) 

FIG. 3. — Error in approximating the cooling (blue lines), heating (red 
lines), and the maximum of the two (black lines) functions with several pho- 
toionization rates. Solid lines trace the median error (50% of cases have the 
error below the solid line), dashed lines trace the 90% error, and dotted line 
trace the 99% error (only 1% of all cases have an error above the dotted line). 
The existence of "catastrophic" errors (a small number of cases with large 
errors) is apparent from this figure. 

MB of memory). 

With these memory requirements, the final table can be 
used even in purely MPI codes on low memory platforms. 

4. TESTING THE COMPLETE APPROXIMATION 

Since we use our sample of Cloudy runs to create the ac- 
tual tables with the cooling and heating functions, we need a 
different data set to test the accuracy of our approximations. 
For this purpose we select 100,000 points from within our pa- 
rameter space (O, sampled uniformly on a logarithmic scale 
(for the metallicity, we randomly choose a value between -3 
and 0.5 in lg(Z)). For each test point, we run Cloudy for our 
8 1 values of the temperature to compute cooling and heating 
functions. This "testing" data set is completely independent 
of the data set used to create the tables. 

We show in Figure [3] the error distribution for our primary 
implementation described above. While errors of both func- 
tions can be evaluated, in practice only the difference between 
the heating and cooling functions matters (EquationQ]) - i.e., 
if one of the two functions is much larger than the other, a 
bigger error can be tolerated for the smaller of the two. Thus, 
we also show in Figure [3] the error of the largest of the two 
functions at each temperature - that error is a more suitable 
measure of the actual error of our approximation. 

Two features of Fig.[3]are important to note. First, the me- 
dian errors for both cooling and heating functions are modest, 
less than 10%. This is very good news indeed, as it shows 
that the whole diversity of cooling and heating functions can 
be parametrized economically, albeit approximately. Second, 
unfortunately, is that the error distribution is not Gaussian, but 
rather exhibits a long tail toward large, or "catastrophic", er- 
rors. For example, in 0.1% of all cases that we tested, the 
error of our approximation reaches a factor of 2. 

This property of our approximation is further illustrated in 
Figure|4] where we show the cumulative error distribution for 
all temperature values separately for the cooling function, the 
heating function, and the maximum of the two (the most ap- 
propriate measure of the actual error). For one case in a mil- 
lion, our approximation reaches an error of about a factor of 
6. Of course, the specific shapes of the distributions shown in 
Fig. [3] and |4] are only applicable to our adopted uniform sam- 
pling. For a specific numerical simulation, the probability of 
errors of a particular magnitude will depend on the simula- 
tion details (such as temperature, density, metallicity, and the 
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FIG. 4. — Cumulative error distributions for all values of gas temperatures 
for the heating (red line), cooling (blue line), and the maximum of the two 
(black line) functions for our approximation. 




T (K) 

FIG. 5. — Cooling (blue lines) and heating (red lines) functions for our test 
models that maximize the error in the cooling function (top panel) and the 
heating function (bottom panel). Approximate functions extracted from our 
table are shown as dashed lines, while actual calculation from Cloudy are 
shown with solid lines. 

radiation field PDFs) and cannot be predicted a priori. 

In order to further illustrate the appearance of a "catas- 
trophic error", we show in Figure [5] the actual cooling and 
heating functions that maximize the errors (out of 100,000 
samples of 81 temperature values each) in both the heating 
and cooling functions. These errors occur in different regions 
of the parameter space (i.e. the worst-case heating function 
is not the heating function that corresponds to the worst-case 
cooling function, and vice versa). In most cases that we were 
able to examine manually, the largest errors occur at either 
very low or very high temperatures, very far from the thermal 
equilibrium values. 

More good news is that the large errors typically occur 
within a narrow range of temperatures, like in Fig. [5] Hence, 
as soon as the gas temperature changes in the simulation code, 
the error in the cooling and heating functions is likely to fall 
appreciably. For example, in the bottom panel of Fig. [5] the 
error in the heating function is a factor of 6 for T = 10 K. The 
heating function is very large, so the gas in this state is going 
to get heated to higher temperature rapidly, trying to reach 
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FIG. 6. — Cumulative error distributions for all values of gas temperatures 
for the maximum of the cooling and heating functions for other spectral 
shapes that are built-in into Cloudy: Milky Way ISM (Cloudy command 
"table ISM", dotted line), typical AGN spectrum (Cloudy command "ta- 
ble AGN", short-dashed line), Haardt-Madau 2005 cosmic background for 
z < 4 (Cloudy command "table HM05", long-dashed line), Haardt-Madau 
2005 cosmic background for z > 4 (Cloudy command "table HM05", dot- 
short-dashed line), continuum from the Crab nebula (Cloudy command "ta- 
ble Crab", dot-long-dashed line), and continuum from a cooling flow (Cloudy 
command "table Cooling Flow", short-dash-long-dashed line). We also show 
with the solid line the black line from Fig.f4]for comparison. 

the thermal equilibrium at T « 2 x 10 6 K. As soon as the gas 
temperature increases to T ~ 100 K, the error in our approx- 
imation drops to less than a factor of 2 (for the fixed values 
of the normalized photoionization rates Qj), and will remain 
within that range until the thermal equilibrium is reached. 

While our adopted spectral shape (0 is sufficiently reason- 
able, it does not represent many situations of astrophysical in- 
terest. In order to test how our approximation fares for other 
radiation fields, we apply it to several spectral shapes that 
are built-in into Cloudy. Specifically, for each of the built- 
in spectral shapes, we randomly choose a value of the density 
within our full range -6 < lg(n&/ cm -3 ) < 6, randomly scale 
the radiation field by a factor between 10~ 3 and 10 3 (except 
for the Haardt-Madau 2005 background, for which we ran- 
domly choose a value of cosmic redshift between and 6), 
use Cloudy to compute the actual cooling and heating func- 
tions for that set of parameters, and then compare Cloudy re- 
sults with the cooling and heating functions returned by our 
approximation when we use the actual photoionization rates 
that Cloudy reports for the chosen spectral shape. 

The cumulative error distribution for several built-in spec- 
tral shapes is shown in Figure [6] In all cases except one the 
maximum error of our approximation does not exceed a factor 
of 3, and the error distribution falls off sharply at large errors 
- in fact, more sharply than the error distribution for our cho- 
sen spectral shape from Equation (O; this is not an artifact 
of incomplete sampling - we run enough Cloudy models for 
each spectral shape to sample the error distribution all the way 
down to 10" 6 . 

The only case where our approximation fares worse, pro- 
ducing a significant fraction (10~ 3 ) errors as large as a fac- 
tor of 5, is the case of the Haardt-Madau 2005 cosmic back- 
ground for z > 4. This is rather disappointing, as cosmolog- 
ical simulations are intended to be the primary user of our 
approximation, but we have not succeeded in finding an ap- 
proximation that has a noticeably lower error for that case. 

5. CONCLUSIONS 
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Our main result is that one can approximately represent the 
most general cooling and heating functions for gas in ioniza- 
tion equilibrium as 

{r 1 A}(r,/i ll z,; 1( )BVf^) {r,A},.(r, ,« ft ), (io) 
,=o V z ©/ 

where rj are given in Equation ©. This approximation is 
rather accurate on average, but suffers from "catastrophic" er- 
rors - in 10~ 6 of all cases the approximate cooling or heating 
function may deviate from the exact calculation by up to a 
factor of 6. Thus, our approximation is not suitable for all 
applications. 

Equation ( fTOb does capture the qualitative dependence of 
the cooling and heating functions on the incident radiation 
field. To illustrate this, we show in the appendix three ex- 
amples where the cooling and heating functions are signifi- 
cantly modified by the incident radiation field. The last exam- 
ple - the quasar irradiating its own galactic halo OA. 31 ) - not 
only shows a large effect the radiation field can have on the 



cooling/heating functions, but actually presents an alternative 
feedback mechanism: the central black hole suppresses the 
gas accretion from the halo without any additional mechani- 
cal or thermal feedback. 

Our data table and the reader code for it are available online 
at |http : / /astro . uchicago . edu/~gnedin| 

We are grateful to Andrey Kravtsov for enlightening dis- 
cussions and constructive criticism. Comments by the anony- 
mous referee helped us to realize the inadequacy of the origi- 
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part by the DOE at Fermilab, by the NSF grant AST-0908063, 
and by the NASA grant NNX-09AJ54G. The calculations 
used in this work have been performed on the Joint Fermilab - 
KICP Supercomputing Cluster, supported by grants from Fer- 
milab, Kavli Institute for Cosmological Physics, and the Uni- 
versity of Chicago. We acknowledge the use of code Cloudy 
dFerland et al.l 1998b as the primary research tool of this study. 
This work made extensive use of the NASA Astrophysics 
Data System and arXiv . org preprint server. 



APPENDIX 

SOME EXAMPLES OF COOLING AND HEATING FUNCTIONS IN ISM AND IGM 

In this section we present a few examples where the incident radiation field significantly affects the cooling and heating rates 
in the gas. These examples are not real physical models, but are simple demonstrations that the dependence that we explore in 
this paper actually matters. 

The examples presented here are not exhaustive, of course; one can imagine many other similar situations. Their purpose is 
to illustrate the numerous possible feedback effects in interstellar and intergalactic environments that arise when we take into 
account the effects of the gas metallicity and the incident radiation field on the cooling and heating rates. These effects can be 
studied, even if only approximately, with the approximations presented in this paper. 

Galactic Halo Near a Quasar 




T (K) T (K) 



FIG. 7. — Left: Cooling (blue lines) and heating (red lines) functions for a Z = Zq, nj, = 340 X «/, galactic halo at the specified distances from a quasar with an 
ionizing luminosity 10 13 Lq. The black solid line shows the pure CIE "standard cooling function". Right: Cooling (blue lines) and heating (red lines) functions 
for a Z = Zq, nt = 1 cm -3 HII region around an O star. The black solid line shows the pure CIE "standard cooling function". 

In the left panel of Figure [7] we show cooling and heating functions for a typical galactic halo at z = {rib = 340 x fib = 
8.5 x 10~ 5 cm" 3 ) surrounding a bright quasar with an ionizing luminosity of 10 13 L Q (roughly corresponding to a 1O 9 M black 
hole). We assume solar metallicity, a quasar spectrum of J v oc v~ 2 , and the lHaardt & Madaul (120011) background. 

Some interesting consequences may arise from the radiation-field-dependence of the cooling and heating functions. For ex- 
ample, gas in the halo within 1 Mpc of this quasar will not be able to cool and condense into the disk if its virial temperature is 
below about 10 5 K. 

HII Region Around an O Star 

In the right panel of Figure [7] we show the cooling and heating functions for a solar metallicity cloud with density rib = 1 cm -3 
surrounding an O star with bolometric luminosity L = 30,000Lq. For the stellar spectrum, we assume a black-body with T = 
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30,000 K. The distances we consider are well within the star's Stromgren radius (~ 30 pc), so we may safely assume that the 
radial dependence of the starlight is 1/r 2 (no depletion due to recombinations). If, instead of a single star, we consider a cluster 
ofNO stars, our result will still hold if we simply rescale the distance axis by N l l 2 . 
Close enough to the star, the equilibrium temperature of the HII region can be substantially higher than the canonical 10 4 K. 

Quasar Irradiating its own Halo 



l i q 
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FIG. 8. — Cooling (blue lines) and heating (red lines) curves for a halo with a z = 3 isothermal density profile surrounding a quasar with 10 13 Lq in ionizing 
radiation. In the left panel, the metallicity is fixed at O.5Z0, in which case the cooling and heating functions become distance-independent. In the right panel, the 
metallicity has a mild outward gradient, Z oc r~ l l 2 . The black solid curves show the pure CIE "standard cooling function". 



In Figure [8] we show cooling and heating functions in a gaseous halo at z = 3 (virial density «^ = 200 x rib = 3.2 x 10 
irradiated by a ~ 10 9 M Q central black hole (10 13 L Q in ionizing radiation). The density profile of the cloud is taken as 

o„ ,^-3 - 3 /100kpc 
n b = 3.2 x 10 cm i I — 

and the metallicity is taken either to be constant 0.5Z Q (leading to distance-independent heating and cooling functions) or to have 
a mild outward gradient, 



Z = 0.5Z, 



'0 



lkp C X 1/2 



In both cases, the quasar is capable of maintaining the heating rate in excess of the cooling rate for T < 10 5 K. It is therefore 
possible to prevent cooling in the halo - and hence, accretion of fresh gas onto the galactic disk and the black hole - without 
any need for a mechani cal or non-radiative thermal feedback mechanism. This result is, of course, not new (Sa/onov et al.ll2.Q05l: 
ICiotti & O striker 20Q3), here we simply use it as an illustration to the importance of properly accounting for the radiation field 
dependence of the cooling and heating functions. 



REFERENCES 



Benjamin, R. A., Benson, B. A., & Cox, D. P. 2001, ApJ, 554, L225 
Benson, A. J., Lacey, C. G., Baugh, C. M., Cole, S., & Frenk, C. S. 2002, 

MNRAS, 333, 156 
Boehringer, H. & Hensler, G. 1989, A&A, 215, 147 
Ciotti, L. & Ostriker, J. P. 2007, ApJ, 665, 1038 
Cox, D. P. & Tucker, W. H. 1969, ApJ, 157, 1 157 
Draine, B. T. 1978, ApJS, 36, 595 

Ferland, G. J., Korista, K. T, Verner, D. A., Ferguson, J. W., Kingdon, J. B., 

& Verner, E. M. 1998, PASP, 110, 761 
Gaetz, T. J. & Salpeter, E. E. 1983, ApJS, 52, 155 
Gnat, O. & Sternberg, A. 2007, ApJS, 168, 213 

Haardt. F. & Madau, P. 2001, in Clusters of Galaxies and the High Redshift 

Universe Observed in X-rays, ed. D. M. Neumann & J. T. V. Tran 
Kravtsov, A. V. 2003, ApJ, 590, LI 
Landi, E. & Landini, M. 1999, A&A, 347, 401 



Leitherer, C, Schaerer, D., Goldader, J. D., Delgado, R. M. G., Robert, C, 
Kune, D. F, de Mello, D. E, Devost, D., & Heckman, T. M. 1999, ApJS, 
123, 3 

Mathis, J. S., Mezger, P. G.. & Panagia, N. 1983, A&A, 128, 212 
Raymond, J. C. Cox, D. P., & Smith. B. W. 1976, ApJ, 204, 290 
Ricotti, M., Gnedin, N. Y., & Shull, J. M. 2002, ApJ, 575, 33 
Santoro, F. & Shull, J. M. 2006, ApJ, 643, 26 

Sazonov, S. Y, Ostriker, J. P., Ciotti, L., & Sunyaev, R. A. 2005, MNRAS, 
358, 168 

Shull, J. M. & van Steenberg, M. 1982, ApJS, 48, 95 
Smith, B., Sigurdsson, S., & Abel, T. 2008, MNRAS, 385, 1443 
Sutherland, R. S. & Dopita, M. A. 1993, ApJS, 88, 253 
Vasiliev, E. O. 2011, MNRAS, 414, 3145 

Wiersma, R. P. C, Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99 



